scRNA-seq
## load required packages
library(Seurat)
library(cowplot)
library(dplyr)
library(ggplot2)
library(DT)
library(paletteer)
library(forcats)DEG Analysis overview
Differentially Expressed genes between TRT and CTL
Number of Cells
res = "sample"
obj.tmp@meta.data %>% ggplot(aes(!!sym(res), fill=!!sym(res))) +
geom_bar(alpha=0.7, color="grey5", size=0.1) +
geom_text(stat="count", aes(label= ..count..), vjust=-0.5, size=3) +
orig.ident_fils +
xlab("") +
theme_classic() +
theme(legend.title = element_blank(),
axis.text.x = element_text(angle = 45, hjust=1)) 
Analysis
Volcano plot
# Define funtion
#id1 = "treatment"
#id2 = "control"
#logfc = 0
deg.two.groups = function(obj.srt = obj.srt, treatment, control,logfc=0){
Idents(obj.srt) = 'sample'
markers <-FindMarkers(
obj.srt,
logfc.threshold = logfc,
ident.1 = treatment,
ident.2 = control, slot= "data")
}deg = deg.two.groups(obj.srt = obj.tmp,
treatment = "TRT",
control = "CTL")
t= paste0(paste0("TRT"), " / ",paste0("CTL"))deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) +
geom_point(size=0.1) +
geom_vline(xintercept = 0, size=0.1) +
theme_bw() + ggtitle(t)
Volcano plot with DEG information
deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
"no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) +
geom_point(size=0.5, shape=19, alpha=0.7) +
geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
geom_hline(yintercept = -log10(0.05), size=0.1) +
scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
theme_bw() +
annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up),
hjust = 1.1, vjust = 2, size = 5, color = "red") +
annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn),
hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
theme_bw() + ggtitle(t)
significance : p value <= 0.05
log2FC >= 1.2
DEG table
# GSEA related functions
library(clusterProfiler)
perform_GSEA <- function(res, ref, pvalueCutoff = 1) {
ranking <- function(res) {
# Check the name of log2fc related
if ("avg_log2FC" %in% names(res)) {
df <- res$avg_log2FC
} else if ("log2FoldChange" %in% names(res)) {
df <- res$log2FoldChange
} else {
stop("Neither avg_log2FC nor log2FoldChange columns found in the data frame.")
}
names(df) <- rownames(res)
df <- sort(df, decreasing = TRUE)
return(df)
}
ranked.res <- ranking(res)
x <- clusterProfiler::GSEA(geneList = ranked.res,
TERM2GENE = ref,
pvalueCutoff = pvalueCutoff,
pAdjustMethod = "BH",
verbose = TRUE,
seed = TRUE)
result <- x@result %>% arrange(desc(NES))
result <- result[, c('NES', 'pvalue', 'p.adjust', 'core_enrichment', 'ID')]
return(result)
}
# GSEA Plot
gsea_nes_plot <- function(gsea.res, title, color="pvalue") {
gsea.res = gsea.res %>% mutate(sig=ifelse(pvalue <= 0.05, "p value <= 0.05", "p value > 0.05"))
# basic plot
p <- gsea.res %>%
ggplot(aes(reorder(ID, NES), NES)) +
geom_col(aes(fill=!!sym(color)), color="grey5", size=0.15, alpha=0.8) +
coord_flip() +
labs(x="Pathway", y="Normalized Enrichment Score", title="GSEA") +
theme_classic() +
theme(axis.text.x = element_text(size=5, face = 'bold'),
axis.text.y = element_text(size=6, face = 'bold'),
axis.title = element_text(size=10)) +
ggtitle(title)
# color by color input type
if (color == "pvalue") {
p <- p + scale_fill_gradient(low = 'orangered', high = '#E5E7E9')
} else if (color == "sig") {
p <- p + scale_fill_manual(values = c("orangered", "#E5E7E9"))
}
return(p)
}
